Getting started:

## install or load package
# install.packages("oceanmap")
library("oceanmap")

## Package overview:
?oceanmap 
help(package="oceanmap") ## list of functions

Introduction

Visualizing data is a crucial step in analyzing and exploring data. During the last two decades the statistical programming language R has become a major tool for data analyses and visualization across different fields of science. However, creating figures ready for scientific publication can be a tricky and time consuming task. The oceanmap package provides some helpful functions to facilitate and optimize the visualization of geographic and oceanographic data, such as satellite and bathymetric data sets. Its major functions are

• plotmap() and ggplotmap() to plot land masks
• set.colorbar(), to define the colorbar (position, ticks, colormap)
• v(), an optimized tool for creating image or contour plots from different data formats (matrix, array, raster, netcdf, binary)
• figure(), to facilitate resizing and saving of figures in diverse formats (jpeg, png, eps, pdf and eps)

These functions were written in a way that they do not require a large amount of their numerous arguments to be specified but still return nice plots. Each of these functions will be subsequently discussed. A closely related problem represents the extraction of spatial data at defined positions/areas. Some examples on this will be discussed in an extra section.

The plotmap-function

Often it is required to add a land mask as geographical reference to an image- or scatter- plot. To add a landmask it is required to provide information about the region extent. As for the entire oceanmap-package, the idea of the plotmap-function is to reduce the user effort when creating figures. To allow a high level of user flexibility plotmap can produce a may by 5 different ways. By providing the required region extent as:

1. geographical coordinates (longitude and latitude),
2. a raster-object
3. an extent-object
4. a region-keyword
5. or add a land mask to an existing plot

Many additional arguments can passed by the user. The idea of a region-keyword (option 4) is to use a short character string instead of requiring axes limits. This keyword feature is also implemented in the v-function and permits here the rapid visualization of subregions of interest (e.g. from a raster object). Region keywords can be set up by the add.region-function. The latter function is interactive and also defines the default colorbar placement that is further used by the v-function.

lon <- c(-6, 37)
lat <- c(30, 46)
plotmap(lon=lon, lat=lat, main="Mediterranean Sea")
plotmap(xlim=lon, ylim=lat, main="Mediterranean Sea")

library('raster')
ext <- extent(lon, lat)
plotmap(ext, main="Mediterranean Sea") # extent-object

r <- raster(ext)
plotmap(r, main="Mediterranean Sea") # raster-object

# run regions() to check pre-installed region labels
# (we will introduce later how to define new region definitions)
plotmap('med4', main="Mediterranean Sea") # region label for the Mediterranean

plot(3.7008, 43.4079, xlim=lon, ylim=lat)
plotmap(add=T)
points(3.7008, 43.4079, pch=19)

par(mfrow=c(2, 1))
plotmap('medw4', main="Western Mediterranean Sea",col.bg="darkblue")
plotmap('medw4', main="Western Mediterranean Sea", bwd=3, border='grey', 
        grid=FALSE)

do.save <- FALSE ## open a plotting window
figure("Gulf_of_Lions_extended", do.save=do.save, width=5, height=5,type="pdf")
plotmap("lion",col.bg='darkblue' ,grid=FALSE)
close_fig(do.save)

width <- dev.size() [1]
height <- dev.size() [2]

do.save <- TRUE ## do NOT open a plotting window, but save figure internally

figure("Gulf_of_Lions_extended", do.save=do.save, width=width, height=height, type="pdf")
plotmap("lion",col.bg='darkblue' ,grid=FALSE)
close_fig(do.save)

par(mfrow=c(2,1))
plotmap(lon=c(80, -120), lat=c(-50, 10), main="map from East to West")
plotmap(lon=c(-120, 80), lat=c(-50, 10), main="map from West to East")

The set.colorbar-function

Setting up the colorbar for an image plot in R is particularly tricky. Common plot functions, such as image.plot of the fields-package set the colorbar to a defined margin, but lack manipulation procedures. Setting up a colorbar by hand requires a lot of data input and is hence time consuming. In particular four different arguments are needed:

• the limits of the colorbar, defining its position, height and width
• the axis where to set tick labels (bottom, left, top, right)
• the colormap and its gradient (horizontal or vertical)
• the labels of the colorbar

The set.colorbar-function of the oceanmap-package comes with 3 ways to define the colorbar position:

1. defining its spatial extent (arguments cbx, cby)

2. defining argument cbpos with a single letter (“b”, “I”, “t”, "r) that indicates the position of the colorbar (bottom, left, top, right)

3. set the position of colorbar manually with the mouse cursor (run set.colorbar() without defining cbx, cby or cbpos).

n either case, the spatial extent of the thus defined colorbar position and extent will be returned (as cbx and cby) that can be reused in later function calls. This feature is further included in the add.region- and v-functions to set up the (default) colorbar placement of regions.

Note that the colormap and ticks of the colorbar can be defined through additional arguments of set.colorbar()

par(mar=c(8,8,8,8))
plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1))
set.colorbar(cbx=c(0, 1), cby=c(-.3, -.4)) # bottom
set.colorbar(cby=c(0, 1), cbx=c(-.4, -.3)) # left
set.colorbar(cbx=c(0, 1), cby=c(1.2, 1.3)) # top
set.colorbar(cby=c(0, 1), cbx=c(1.2, 1.3)) # right

par(mar=c(8,8,8,8))
plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1))
set.colorbar(cbpos='b') # bottom
set.colorbar(cbpos='l') # left
set.colorbar(cbpos='t') # top
set.colorbar(cbpos='r') # right

par(mar=c(8,8,8,8))
plot(0.5,0.5,xlim=c(0,1),ylim=c(0,1))
cb <- set.colorbar() # interactive
plot (0.5,0.5,xlim=c(0,1),ylim=c(0,1))
set.colorbar(cbx=cb$cbx, cby=cb$cby) # reuse stored colorbar

The v-function: Creating image-plots of spatial data

The v-function combines the advantages of the plotmap and set.colorbar to visualize 2D (oceanographic) data. To facilitate plotting, it is further capable of handling different input formats. Valid formats are:

• matrix and array objects
• objects of class ‘Raster’ (‘RasterLayer’, ‘RasterStack’ or ‘RasterBrick’),
• netcdf-files ‘nedf4’-objects,
• ‘.gz’-files (Hervé Demarq)

Note that for plotting the input data will be internally transformed into a Raster-object, no matter which of the valid input format was selected. Examples to visualize matrix and array objects will not be discussed here but are given in ?matrix2raster(). The other input data formats are discussed in the following sections.

setwd(system.file("test_files", package="oceanmap"))
load("medw4_modis_sst2_4km_1d_20020705_20020705.r2010.0.qual0.Rdata")
dat <- raster::crop(dat,extent(c(0,10,40,44))) ## crop data, xlim/ylim not yet implemented in v()
print(dat)
## class      : RasterLayer 
## dimensions : 96, 240, 23040  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : 0, 10, 40, 44  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=WGS84 
## source     : memory
## names      : layer 
## values     : 16.8, 25.8  (min, max)

v(dat, main="Raster-object", cbpos='r',subplot=TRUE)

setwd(system.file("test_files", package="oceanmap"))
ncfiles <- Sys.glob('*.nc') # load list of sample-'.nc'-files
print(ncfiles)

## [1] "herring_lavae.nc" "topo.v02.nc"

library('ncdf4')
ncdf <- nc_open(ncfiles[1])
print(ncdf)

## File herring_lavae.nc (NC_FORMAT_CLASSIC):
## 
##      1 variables (excluding dimension variables):
##         float Conc[lon,lat,time]   
##             _FillValue: -1
## 
##      3 dimensions:
##         lon  Size:352
##         lat  Size:310
##         time  Size:4   *** is unlimited ***
##             units: seconds since 2006-03-01 00:00:00

par(mar=c(5,5,5,5))
v(obj = ncdf, cbpos="r",subplot=TRUE)

setwd(system.file("test_files", package="oceanmap"))
nc <- nc2raster(ncfiles[1])

## File herring_lavae.nc (NC_FORMAT_CLASSIC):
## 
##      1 variables (excluding dimension variables):
##         float Conc[lon,lat,time]   
##             _FillValue: -1
## 
##      3 dimensions:
##         lon  Size:352
##         lat  Size:310
##         time  Size:4   *** is unlimited ***
##             units: seconds since 2006-03-01 00:00:00
## 
## converting array to RasterStack object

setwd(system.file("test_files", package="oceanmap"))
## attention you may have to increase margin spacings for correct colorbar text rendering
par(mar=c(5,5,5,5)) 
v(nc, cbpos="r",subplot=TRUE) # plot RasterStack object

setwd(system.file("test_files", package="oceanmap"))
par(mfrow=c(2,1),mar=c(5,5,5,5))
v(ncfiles[1], cbpos="r", subplot = TRUE)
v(ncfiles[1], cbpos="r", replace.na=TRUE, subplot = TRUE)

setwd(system.file("test_files", package="oceanmap"))
par(mfrow=c(2,2),mar=c(5,5,5,5))
v(ncfiles[1], t=1:4, cbpos="r", subplot = TRUE)

lon <- c(9, 31)
lat <- c(53.5, 66)
get.bathy(lon=lon, lat=lat, main="Baltic Sea", cbpos='r', subplot = TRUE)

setwd(system.file("test_files", package="oceanmap"))
par(mfrow=c(2,1))
bathy <- get.bathy("medw4", terrain=T, res=10, keep=T, visualize=T, subplot = TRUE, grid=F)
# load('bathy_medw4_res.3.dat',verbose = T); bathy <- h
v(bathy, param="bathy", subplot = TRUE, terrain=F, levels=c(200 ,2000)) # show contours

par(mfrow=c(2,1))
h <- get.bathy("lion",terrain=FALSE, res=10, visualize=TRUE, 
               subplot=TRUE, v_image = FALSE, levels=c(200,2000) )

## use v-function for same plot but on subregion:
v(h,v_area = "survey", param="bathy", subplot=TRUE, 
  v_image = FALSE, levels=c(200,2000))

v() and gz-files:

owd <- getwd()
setwd(system.file("test_files", package="oceanmap"))
gz.files <- Sys.glob('*.gz')
# v(gz.files[2],subplot=F) ## plot content of gz-file

# obj <- readbin(gz.files[2],area='lion')
# par (mfrow=c(1,2))
# v(obj ,param="sst",subplot = TRUE)
# v(obj,param="Temp",subplot = TRUE) ## note unset "pal" (colormap) for unkown "param"-values!

available colormaps

oceanmap comes with a set of pre-installed colormaps. However, the user can apply own colormaps in v-function calls.

data('cmap') # load color maps data
names(cmap) # list available color maps

##  [1] "ano"            "bathy"          "blue"           "chla"          
##  [5] "green"          "haxby"          "jet"            "rainbow"       
##  [9] "red"            "sst"            "haxbyrev"       "orange"        
## [13] "year.jet"       "dark.blue"      "jetrev"         "light.jet"     
## [17] "white2black"    "black2white"    "white2darkgrey"

setwd(system.file("test_files", package="oceanmap"))
gz.files <- Sys.glob('*.gz')
# figure(width=15,height=15)
par(mfrow=c(4,5))
# for(n in names(cmap)) v(gz.files[2], v_area='lion', subplot=TRUE, pal=n, adaptive.vals=TRUE, main=n)
## define new color maps from blue to red to white:
n <- colorRampPalette(c('blue','red','white')) (100)
# v(gz.files[2], v_area='lion', subplot=TRUE, pal=n, adaptive.vals=TRUE, main="own colormap")

available oceanographic parameters

The standard param values required in v() to define the title and colormap of the colorbar are stored in the parameter definitions dataset. Here some examples of pre-installed definitions.

## Example 6: available parameters
data(parameter_definitions)
names (parameter_definitions)

##  [1] "param"           "a"               "b"               "c"              
##  [5] "log"             "name1"           "unit"            "pal1"           
##  [9] "minv"            "maxv"            "min"             "max"            
## [13] "invalid_data_dc" "coast_dc"        "land_dc"         "no_data_dc"

# ?parameter_ definitions
setwd(system.file("test_files", package="oceanmap"))

# figure (width=12, height=6.2)
par(mfrow=c(2,3))
# v('*sst2*707*' ,v_area="medw4",main="sst", subplot=TRUE)
# v('*chla*531*' ,v_area="medw4" ,main="chla", subplot=TRUE)
# v('*chlagrad*' ,v_area="medw4" ,main="chlagrad", subplot=TRUE)
# v('?*p100*' ,v_area="medw4",main="p100 (oceanic fronts)",subplot=TRUE)
# v('*xsla*',v_area="medw4",main="sla",subplot=TRUE)
# h <- get.bathy("medw4",visualize=TRUE,terrain=F,res=4, subplot=TRUE, main="bathy")

available region keywords & definitions

Researchers often need to analyse spatial data of a fixed study region (e.g. the Western Mediterranean Sea). To facilitate plotting of such data, oceanmap comes with a region- keyword feature. The idea of this is to link all information required to produce land mask or image plots to a short keyword (label). Apart from its keyword, each region-definition includes a long version of the label, the spatial extent of the region (longitudes and latitudes), its grid resolution, as well as the default colorbar position and default figure dimensions. Available region-definitions are stored as a data frame in the region_definitions file of the package.

data(region_definitions)
head (region_definitions)
region_definitions$label
# ?region_definitions
# figure (width=15,height=15)
par(mfrow=c(7,4))
for(n in region_definitions$label) plotmap(region = n,main=n)

How to create and delete region-definitions with add.region() - part I

The function add.region() can be used to set up new region-definitions (and -keywords), to store or restore a backup. ATTENTION: Backups are necessary since new region- definitions will be lost when updating eceanmap! In general, new regions can be supplied as a one-row data frame, e.g. by modifying existing region definition entries, or as a second option, through an interactive process. Here an example on how to modify existing definitions (the interactive way will be discussed on the next page):

# data(region_definitions) # load region_definitions
# lion <- region_definitions[region_definitions$label == 'lion',] # selecting Gulf of Lions region

# delete.region("junk") # delete junk region
# data(region_definitions) # reload region_definitions
# region_definitions[,1:9]

How to create new region-definitions with add.region() - part II: the interactive way

Calling add.region() without or insufficient arguments starts an interactive session that prompts the user to complete the new region definition. The single steps of this way listed below:

Please define the keyword of the new region, coded as 'label':
Please define the long name of the new region, coded as 'name':

Please define the northern most latitude
(negative values for the southern hemisphere) of the new region, coded as 'latn':

Please define the southern most latitude
(negative values for the southern hemisphere) of the new region, coded as 'lats':

Please define the western most longitude of the new region
(negative values for the western hemisphere), coded as 'lonw':

Please define the western most longitude of the new region
(negative values for the western hemisphere), coded as 'lone':

Please type 'm' if you want to perform colorbar placement by hand (mouse cursor)
or a letter (b, 1, t, r) referring to a side of the plot (bottom, left, top, right).

If 'm' (manual placement) was selected:
Please select the lower left colorbar-position by the mouse cursor.

Please select the upper right colorbar-position by the mouse cursor.

You can resize the window to appropriate size. Try to avoid white space.
Press <Enter> when done.

Please enter the default grid resolution

Press <Enter> to save the new region definition
or any other key to abort the operation

Extra: Extracting data from spatial Raster-objects:

Extracting spatial data from Raster-objects can be easily done with a bunch of functions implemented in the sp and maptools packages (in particular with the function extract). Here some examples, for points, polygons and circles:

bathy <- get.bathy('lion',subplot=TRUE,res=10) # get bathymetry

## loading bathymetry data at a resolution of 10 minutes
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size

## 
## printing file: lion_bathy

#' extract values
# Points <- locator()
Points <- cbind(lon=c(6,8),lat=c(40,40))
points(Points,pch=19,col=c('black','red'))
extract(bathy, Points) # extract points

## [1] 2829 1379

#' polygon
library('rgeos')
cds1 <- rbind(c(3,42), c(3,43), c(7,43), c(3,42)) # polygon need to be closed (first point = last point)
poly <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1)))
plot(poly[1],add=T,lwd=2,border="blue")
extract(bathy, poly[1])

## [[1]]
##  [1]  NaN   43   72   91  104  107  178  127  194  367  560  423  992 1322 1708
## [16] 1905 2024 2158 2223 1949 2101 2303  NaN   48   83   95  211  257  480  286
## [31]  532  809  976  982 1495 1739 1843 2032 2267 2374  NaN   51   89  132  309
## [46]  624  702  901 1369 1254 1349 1412 1732 1847  NaN  NaN   98  228  719  903
## [61] 1091 1174 1694 1754  NaN    1   47  126  764 1314  NaN  NaN

#
poly2 <- gBuffer(poly[1], width = 10/60) # increase polygon by 10 nm
plot(poly2,add=T,lwd=2,border="blue",lty='dotted')

# source("SpatialCircle.r")
library(maptools) 
circ <- SpatialCircle(Points[1,1],Points[1,2],r = 1)
circ_poly <- PolySet2SpatialPolygons(SpatialLines2PolySet(circ)) # convert circle from spatialline to spatialpolygon
plot(circ_poly, add=T)

extract(bathy, circ_poly)

## [[1]]
##   [1] 2692 2696 2700 2707 2719 2714 2721 2719 2727 2736 2753 2765 2745 2748 2748
##  [16] 2749 2748 2753 2760 2781 2795 2803 2776 2774 2772 2776 2778 2779 2782 2805
##  [31] 2821 2829 2752 2764 2792 2811 2823 2827 2826 2826 2769 2795 2828 2643 2688
##  [46] 2709 2746 2839 2877 2858 2842 2842 2769 2791 2831 2849 2714 2785 2832 2865
##  [61] 2810 2794 2829 2888 2738 2700 2735 2812 2787 2844 2812 2825 2826 2781 2777
##  [76] 2847 2811 2847 2856 2857 2634 2788 2812 2862 2827 2824 2857 2824 2733 2689
##  [91] 2599 2735 2808 2847 2852 2850 2865 2855 2828 2759 2755 2787 2892 2868 2857
## [106] 2919 2918 2849 2909 2861 2828 2852